function wint=interpolation(ewsb, index,gammarfun);
global P beta n m S state agent gammagrid Y scalefactor lowerbound rho table elementsize

lambda(1)=1/(1+gammarfun(1)+gammarfun(2));
lambda(2)=gammarfun(1)/(1+gammarfun(1)+gammarfun(2));
lambda(3)=gammarfun(2)/(1+gammarfun(1)+gammarfun(2));
lambda=lambda(1,1:3)+eps;
if lambda(1:3)>=lowerbound & isreal(lambda);
node=fix((lambda-lowerbound)/elementsize+1);
jj=index(2);

for g=1:3
   if node(g)==0
       node(g)=1;
   elseif node(g)==m
       node(g)=m-1;
   end
end

x(1)=[(node(1)-1)*elementsize+lowerbound];
x(2)=x(1)+elementsize;
y(1)=[(node(2)-1)*elementsize+lowerbound];
y(2)=y(1)+elementsize;
z=ones(2,2,agent)*-9999;
for r1=0:1
    for r2=0:1
    if table(node(1)+r1,node(2)+r2)~=0
z(r1+1,r2+1,:)=ewsb(table(node(1)+r1,node(2)+r2),jj,:);
    end
    end
end
if table(node(1),node(2))==0
z(1,1,:)=z(2,1,:);
end
if table(node(1)+1,node(2)+1)==0
z(2,2,:)=z(1,2,:);
end
for g=1:agent


fR1=(lambda(1)-x(2))/(x(1)-x(2))*z(1,1,g) + (lambda(1)-x(1))/(x(2)-x(1))*z(2,1,g);
fR2=(lambda(1)-x(2))/(x(1)-x(2))*z(1,2,g)+ (lambda(1)-x(1))/(x(2)-x(1))*z(2,2,g);
wint(g)=(lambda(2)-y(2))/(y(1)-y(2))*fR1 + (lambda(2)-y(1))/(y(2)-y(1))*fR2; 

end


else
    wint(1:agent)=-9999;
end
